4. Clustering and classification

In this chapter I study “Boston” dataset from the MASS package. The chapter starts with an introduction to the dataset while the focus is on clustering and classification in which the idea is to figure out clusters of data points that are in some sense closer to each other than some other data points. After finding some clusters, new observations can be classified to these clusters. As discussed later, linear discriminant analysis (LDA) is one way to find and separate clusters from each other: LDA finds the (linear) combination of the variables that separate the target variable classes.

Part of this exercise includes data wrangling for next chapter, which is not commented here. For the data wrangling, see the github repository.

Analysis

The chapter is divided to subsections (Part 2 - Part 7 + bonus). Each subsection has a brief introduction to inform about its content.

I start with some house keeping and calling the packages used later in the exercise.

rm(list=ls()) 
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(GGally)
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa
library(tidyr)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(corrplot)
## corrplot 0.84 loaded

Part 2: Importing data + brief overview

First thing here is to load the data from MASS-package, after which I print out the names of the variables in the data and describe the data set briefly. The “Boston” dataset is described here.

data("Boston")
str(Boston)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
summary(Boston)
##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

The dataset contains information about the housing values in suburbs in Boston with 506 observations and 14 variables. The variables in our dataset are:

  • crim: per capita crime rate by town.

  • zn: proportion of residential land zoned for lots over 25,000 sq.ft.

  • indus: proportion of non-retail business acres per town.

  • chas: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).

  • nox: nitrogen oxides concentration (parts per 10 million).

  • rm: average number of rooms per dwelling.

  • age: proportion of owner-occupied units built prior to 1940.

  • dis: weighted mean of distances to five Boston employment centres.

  • rad: index of accessibility to radial highways.

  • tax: full-value property-tax rate per $10,000.

  • ptratio: pupil-teacher ratio by town.

  • black: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.

  • lstat: lower status of the population (percent).

  • medv: median value of owner-occupied homes in $1000s

Source: Harrison, D. and Rubinfeld, D.L. (1978) Hedonic prices and the demand for clean air. < em >J. Environ. Economics and Management < b >5, 81–102.

Belsley D.A., Kuh, E. and Welsch, R.E. (1980) < em >Regression Diagnostics. Identifying Influential Data and Sources of Collinearity. New York: Wiley.

Part 3: Graphical overview of the data and summaries of the variables

In this subsection the dataset is introduced in a more detailed manner. The variables and their relationships are studied both graphically and numerically.

pairs(Boston)

summary(Boston)
##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00
cor_matrix<-cor(Boston) %>% round(2)
cor_matrix
##          crim    zn indus  chas   nox    rm   age   dis   rad   tax
## crim     1.00 -0.20  0.41 -0.06  0.42 -0.22  0.35 -0.38  0.63  0.58
## zn      -0.20  1.00 -0.53 -0.04 -0.52  0.31 -0.57  0.66 -0.31 -0.31
## indus    0.41 -0.53  1.00  0.06  0.76 -0.39  0.64 -0.71  0.60  0.72
## chas    -0.06 -0.04  0.06  1.00  0.09  0.09  0.09 -0.10 -0.01 -0.04
## nox      0.42 -0.52  0.76  0.09  1.00 -0.30  0.73 -0.77  0.61  0.67
## rm      -0.22  0.31 -0.39  0.09 -0.30  1.00 -0.24  0.21 -0.21 -0.29
## age      0.35 -0.57  0.64  0.09  0.73 -0.24  1.00 -0.75  0.46  0.51
## dis     -0.38  0.66 -0.71 -0.10 -0.77  0.21 -0.75  1.00 -0.49 -0.53
## rad      0.63 -0.31  0.60 -0.01  0.61 -0.21  0.46 -0.49  1.00  0.91
## tax      0.58 -0.31  0.72 -0.04  0.67 -0.29  0.51 -0.53  0.91  1.00
## ptratio  0.29 -0.39  0.38 -0.12  0.19 -0.36  0.26 -0.23  0.46  0.46
## black   -0.39  0.18 -0.36  0.05 -0.38  0.13 -0.27  0.29 -0.44 -0.44
## lstat    0.46 -0.41  0.60 -0.05  0.59 -0.61  0.60 -0.50  0.49  0.54
## medv    -0.39  0.36 -0.48  0.18 -0.43  0.70 -0.38  0.25 -0.38 -0.47
##         ptratio black lstat  medv
## crim       0.29 -0.39  0.46 -0.39
## zn        -0.39  0.18 -0.41  0.36
## indus      0.38 -0.36  0.60 -0.48
## chas      -0.12  0.05 -0.05  0.18
## nox        0.19 -0.38  0.59 -0.43
## rm        -0.36  0.13 -0.61  0.70
## age        0.26 -0.27  0.60 -0.38
## dis       -0.23  0.29 -0.50  0.25
## rad        0.46 -0.44  0.49 -0.38
## tax        0.46 -0.44  0.54 -0.47
## ptratio    1.00 -0.18  0.37 -0.51
## black     -0.18  1.00 -0.37  0.33
## lstat      0.37 -0.37  1.00 -0.74
## medv      -0.51  0.33 -0.74  1.00
corrplot(cor_matrix, method="circle", type="upper", cl.pos="b",tl.pos="d", tl.cex=0.6)

Going through the descriptions of all the distributions and correlations here would be tiresome. To make this easier, pairs() function and summary gives an understanding about the distributions of the variables while correlations between variables are plotted above with corrplot()-function. Correlations always get values between -1 and 1, where 1 refers to a situation in which variables are (perfectly) positively correlated and they move together hand in hand, and -1 refers to a situation where variables move together oppositely: when one decreases the other increases. In the graph above, the bigger the correlation in absolute values, the larger the circle. Red color refers to a negative correlation and blue color refers to a positive correlation. The darker the color, the stronger the correlation. Just to name a few, “tax” and “rad” have a strong positive correlation (0.91) such that the full-value property-tax rate is high in areas with good access to radial highways whereas “nox” and “dis” have a strong negative correlation (-0.77) such that the levels of pollution are high close to the five Boston employment centres.

Part 4: Standardizing the data, creating a categorial variable and train and test sets

Here I standardize the data and study how this affects our variables. I also create a categorical variable of the crime rate and divide the dataset to train and test sets, so that 80% of the data belongs to the train set.

library(dplyr)
#center and standardize variables
boston_scaled <- scale(Boston)
# summaries of the scaled variables
summary(boston_scaled)
##       crim                 zn               indus        
##  Min.   :-0.419367   Min.   :-0.48724   Min.   :-1.5563  
##  1st Qu.:-0.410563   1st Qu.:-0.48724   1st Qu.:-0.8668  
##  Median :-0.390280   Median :-0.48724   Median :-0.2109  
##  Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.007389   3rd Qu.: 0.04872   3rd Qu.: 1.0150  
##  Max.   : 9.924110   Max.   : 3.80047   Max.   : 2.4202  
##       chas              nox                rm               age         
##  Min.   :-0.2723   Min.   :-1.4644   Min.   :-3.8764   Min.   :-2.3331  
##  1st Qu.:-0.2723   1st Qu.:-0.9121   1st Qu.:-0.5681   1st Qu.:-0.8366  
##  Median :-0.2723   Median :-0.1441   Median :-0.1084   Median : 0.3171  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:-0.2723   3rd Qu.: 0.5981   3rd Qu.: 0.4823   3rd Qu.: 0.9059  
##  Max.   : 3.6648   Max.   : 2.7296   Max.   : 3.5515   Max.   : 1.1164  
##       dis               rad               tax             ptratio       
##  Min.   :-1.2658   Min.   :-0.9819   Min.   :-1.3127   Min.   :-2.7047  
##  1st Qu.:-0.8049   1st Qu.:-0.6373   1st Qu.:-0.7668   1st Qu.:-0.4876  
##  Median :-0.2790   Median :-0.5225   Median :-0.4642   Median : 0.2746  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.6617   3rd Qu.: 1.6596   3rd Qu.: 1.5294   3rd Qu.: 0.8058  
##  Max.   : 3.9566   Max.   : 1.6596   Max.   : 1.7964   Max.   : 1.6372  
##      black             lstat              medv        
##  Min.   :-3.9033   Min.   :-1.5296   Min.   :-1.9063  
##  1st Qu.: 0.2049   1st Qu.:-0.7986   1st Qu.:-0.5989  
##  Median : 0.3808   Median :-0.1811   Median :-0.1449  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.4332   3rd Qu.: 0.6024   3rd Qu.: 0.2683  
##  Max.   : 0.4406   Max.   : 3.5453   Max.   : 2.9865
# change the object to data frame
boston_scaled<-as.data.frame(boston_scaled)
# summary of the scaled crime rate
summary(boston_scaled$crim)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.419367 -0.410563 -0.390280  0.000000  0.007389  9.924110
# create a quantile vector of crim and print it
bins <- quantile(boston_scaled$crim)
bins
##           0%          25%          50%          75%         100% 
## -0.419366929 -0.410563278 -0.390280295  0.007389247  9.924109610
# create a categorical variable 'crime' -- I use the bins as breaks
crime <- cut(boston_scaled$crim, breaks = bins, include.lowest = TRUE, labels=c("low", "med_low","med_high","high"))

# remove original crim from the dataset
boston_scaled <- dplyr::select(boston_scaled, -crim)

# add the new categorical value to scaled data
boston_scaled <- data.frame(boston_scaled, crime)

# look at the table of the new factor crime
table(crime)
## crime
##      low  med_low med_high     high 
##      127      126      126      127
# number of rows in the Boston dataset 
n <- nrow(boston_scaled)

# choose randomly 80% of the rows
ind <- sample(n,  size = n * 0.8)

# create train set
train <- boston_scaled[ind,]

# create test set 
test <- boston_scaled[-ind,]

# save the correct classes from test data
correct_classes <- test$crime

# remove the crime variable from test data
test <- dplyr::select(test, -crime)

As shown in the summary of our variables above, scaling the variables changes the values of the observations. In the scaling we subtract the column means from the corresponding columns and divide the difference with standard deviation. As a result, means of the variables become zero and the values vary in a smaller scale. Notice that I had “crim” in original set, but in the LDA-part I use “crime” as new variable. “Crime” is a categorical variable (see code above) that gets values “low”, “med_low”,“med_high” and “high” depending on the values of the original “crim” variable. For these categories, I created bins based on quantiles of the data. Also, the dataset is divided to train and test sets, so that 80% of the data belongs to the train set.

Part 5: Linear discriminant analysis

In this section I fit the linear discriminant analysis on the train set and use the recently created categorical crime rate “crime” as the target variable and all the other variables in the dataset as predictor variables. After that, I draw the LDA (bi)plot.

# linear discriminant analysis
str(train)
## 'data.frame':    404 obs. of  14 variables:
##  $ zn     : num  1.764 0.37 -0.487 -0.487 -0.487 ...
##  $ indus  : num  -0.848 -0.609 -0.437 1.231 -0.437 ...
##  $ chas   : num  -0.272 -0.272 -0.272 -0.272 -0.272 ...
##  $ nox    : num  -1.292 -0.783 -0.144 2.73 -0.144 ...
##  $ rm     : num  0.0432 -0.0635 -0.7937 -1.2534 -0.2685 ...
##  $ age    : num  -0.8162 -1.8571 0.0329 1.1164 1.0063 ...
##  $ dis    : num  1.67267 0.301066 0.000692 -1.046413 -0.016737 ...
##  $ rad    : num  -0.408 -0.752 -0.637 -0.522 -0.637 ...
##  $ tax    : num  -0.6837 -1.0991 -0.6007 -0.0311 -0.6007 ...
##  $ ptratio: num  -0.8571 0.0667 1.1753 -1.7347 1.1753 ...
##  $ black  : num  0.441 0.441 0.375 -0.165 -1.187 ...
##  $ lstat  : num  -0.7076 -0.849 -0.1923 0.0878 1.0764 ...
##  $ medv   : num  -0.0253 0.29 -0.4711 -0.3189 -0.9821 ...
##  $ crime  : Factor w/ 4 levels "low","med_low",..: 1 2 3 3 3 1 1 3 3 2 ...
lda.fit <- lda(crime ~., data = train)
# print the lda.fit object
lda.fit
## Call:
## lda(crime ~ ., data = train)
## 
## Prior probabilities of groups:
##       low   med_low  med_high      high 
## 0.2475248 0.2574257 0.2500000 0.2450495 
## 
## Group means:
##                   zn      indus        chas        nox         rm
## low       1.08420686 -0.9290890 -0.07547406 -0.8882389  0.4518757
## med_low  -0.07124951 -0.2949711 -0.04518867 -0.5452854 -0.1223992
## med_high -0.36667477  0.1677386  0.11748284  0.3838803  0.1730296
## high     -0.48724019  1.0171737 -0.03371693  1.0457030 -0.3730469
##                 age        dis        rad        tax    ptratio
## low      -0.9267125  0.9006450 -0.6947544 -0.7478303 -0.4542995
## med_low  -0.2983417  0.3324553 -0.5434660 -0.4804160 -0.0785039
## med_high  0.4225889 -0.3802532 -0.3780732 -0.2980777 -0.2753546
## high      0.8077492 -0.8427495  1.6375616  1.5136504  0.7801170
##                black       lstat         medv
## low       0.37371838 -0.78328176  0.555087083
## med_low   0.30727155 -0.13301806 -0.003357932
## med_high  0.08799068  0.01494718  0.199359526
## high     -0.75430345  0.82188887 -0.675385876
## 
## Coefficients of linear discriminants:
##                   LD1         LD2           LD3
## zn       0.1200160699  0.82712504 -0.9082491407
## indus   -0.0093357175 -0.23401271  0.1992602384
## chas    -0.0939217255 -0.02031132  0.0848417933
## nox      0.3946105741 -0.61086934 -1.4509535602
## rm      -0.0795084395 -0.15275067 -0.2329164386
## age      0.3048302924 -0.35050291 -0.1405501803
## dis     -0.0716812369 -0.34423163  0.0794969945
## rad      2.9938542244  0.94537879 -0.0001759071
## tax     -0.0007574562 -0.11659670  0.5542646842
## ptratio  0.1242431563  0.06793073 -0.3425059012
## black   -0.1235541110 -0.02952061  0.0686779449
## lstat    0.2061795332 -0.23202443  0.2794142889
## medv     0.1860762414 -0.35753239 -0.2803084799
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.9417 0.0431 0.0152
# the function for lda biplot arrows
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "red", tex = 0.75, choices = c(1,2)){
  heads <- coef(x)
  arrows(x0 = 0, y0 = 0, 
         x1 = myscale * heads[,choices[1]], 
         y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
  text(myscale * heads[,choices], labels = row.names(heads), 
       cex = tex, col=color, pos=3)
}

# target classes as numeric
classes <- as.numeric(train$crime)

# plot the lda results
plot(lda.fit, dimen = 2, col=classes, pch=classes)
lda.arrows(lda.fit, myscale = 2)

Accurate description for the biplot is not given here as it is a topic of next chapter of the IODS-project. However, the biplot suggests that the most influencial linear separators for the clusters are “rad”, “nox” and “zn” as their arrows are the longest in the graph shown above. In short, this means that index of accessibility to radial highways, the nitrogen oxides concentration (parts per 10 million) and the proportion of residential land zoned for lots over 25,000 sq.ft. are the most influencial factors that separate our three classes from each other.

Part 6: Predicting classes with the LDA model

In this part the focus is on predicting the classes of the observations of test data with the LDA-model from previous part. To do this, I have earlier dropped the “crime” variable from the test dataset. Next, I observe how well the LDA model predicts the classes of the test dataset to study the performance of the LDA model.

# predict classes with test data
lda.pred <- predict(lda.fit, newdata = test)

# cross tabulate the results
table(correct = correct_classes, predicted = lda.pred$class)
##           predicted
## correct    low med_low med_high high
##   low       10      16        1    0
##   med_low    1      19        2    0
##   med_high   0      10       15    0
##   high       0       0        0   28

The performance of the LDA model is illustrated in the table above. A quick look shows that most of the observations are on the diagonal elements, which suggests that the model predicts the classes of the test dataset relatively well. To be more precise, the model gets 76 out of 102 predictions correctly (thats approximately 75%) and 26 out of 102 wrong (thats approximately 25%).

It seems that the model has no difficulties categorizing the cases of high crime rate correctly (30/30) while other categories have some troubles in classification. The “low” is the most accurate after “high” with 70% of predictions of real “lows” being categorized correctly while the same numbers for “med_high” and “med_low” are 63% and 60%, respectively. Intuitively, it makes sense that the more extreme cases are easier to categorize correctly than the observations that fall in to the middle categories.

Part 7: Clusters

In this subsection, I reload the original “Boston” dataset, scale it and calculate the distances between the observations, after which run k-means algorithm on the dataset. The goal is to find optimal number of clusters and to visualize the results.

For the distances, I calculate both eucledean and manhattan distances, which are saved in distance martices “dist_eu” and “dist_man”. Notice that the distances are calculated from scaled Boston dataset (“bostonscaled”).

# load MASS and Boston
library(MASS)
data('Boston')
bostonscaled<-scale(Boston)
# euclidean distance matrix
dist_eu <- dist(bostonscaled)

# look at the summary of the distances
summary(dist_eu)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1343  3.4625  4.8241  4.9111  6.1863 14.3970
# manhattan distance matrix
dist_man <- dist(bostonscaled, method="manhattan")

# look at the summary of the distances
summary(dist_man)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2662  8.4832 12.6090 13.5488 17.7568 48.8618

Summaries of the distance matrices show that the mean of eucledean distances is 4.9 (min. 0.13 and max. 14.39) while the mean of manhattan disctances is 13.5 (min. 0.27 and max. 48.7).

Next step here is to run the k-means algorithm to find the optimal number of clusters.

# k-means clustering
set.seed(123)

# determine the number of clusters
k_max <- 10

# calculate the total within sum of squares
twcss <- sapply(1:k_max, function(k){kmeans(Boston, k)$tot.withinss})

# visualize the results
qplot(x = 1:k_max, y = twcss, geom = 'line')

# k-means clustering
km1 <-kmeans(Boston, centers = 2)

# plot the Boston dataset with clusters
pairs(Boston, col = km1$cluster)

The optimal number of clusters is when the total of within cluster sum of squares (total WCSS) drops radically. Studying the qplot suggests that the optimal number could therefore be 2. The qplot show how the total WCSS behaves when the number of cluster changes.

Running the algorithm again and plotting the results gives us the graph of the scaled Boston dataset with 2 clusters (red and black). Interpreting the plot above suggests that compared to the red cluster, the black cluster is categorized by more crime (“crim”), higher levels of pollution (“nox”), closer proximity to railways (“rad”).

Bonus: Finding the most influencial linear separators for the clusters

Here we do the clustering again but with a larger number of clusters (>2) for the scaled Boston dataset. The qplot suggests that if not 2 clusters are used, 3 could be a good choice as well. The line drops quite evenly after 2 so that a larger number of clusters might also be justified. Thus, I set number of clusters to 3 according to qplot.

library(MASS)
data('Boston')
bostonscaled<-scale(Boston)
bostonscaled<-as.data.frame(bostonscaled)


# target classes as numeric
classes <- as.numeric(train$crime)
# k-means clustering
set.seed(123)

# determine the number of clusters
k_max <- 10

# calculate the total within sum of squares
twcss <- sapply(1:k_max, function(k){kmeans(bostonscaled, k)$tot.withinss})
# visualize the results
qplot(x = 1:k_max, y = twcss, geom = 'line')

# k-means clustering
km <-kmeans(bostonscaled, centers = 3)

# plot the Boston dataset with clusters
pairs(bostonscaled, col = km$cluster)

Graph above illustrates the differences between the clusters. Now that the number of clusters is set to 3, I run the LDA again using the clusters as target classes and including all the variables in the Boston dataset to the model.

lda.fit2 <- lda(km$cluster ~., data = bostonscaled)

# the function for lda biplot arrows
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "red", tex = 0.75, choices = c(1,2)){
  heads <- coef(x)
  arrows(x0 = 0, y0 = 0, 
         x1 = myscale * heads[,choices[1]], 
         y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
  text(myscale * heads[,choices], labels = row.names(heads), 
       cex = tex, col=color, pos=3)
}




# plot the lda results
plot(lda.fit2, dimen = 2, col=classes, pch=classes)
lda.arrows(lda.fit2, myscale = 3)

Fitting the LDA model again and studying the biplot suggests that the most influencial linear separators for the clusters are “rad”, “tax” and “age” as their arrows are the longest in the graph shown above. In short, this means that index of accessibility to radial highways,full-value property-tax rate per $10,000 and the proportion of owner-occupied units built prior to 1940 are the most influencial factors that separate our three classes from each other.

Super bonus: matrix product and 3D-plot

In this part, I create a matrix product, which is a projection of the data point and draw a 3D plot where the color is defined by the clusters of the k-means.

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
model_predictors <- dplyr::select(train, -crime)
# check the dimensions
dim(model_predictors)
## [1] 404  13
dim(lda.fit$scaling)
## [1] 13  3
# matrix multiplication
matrix_product <- as.matrix(model_predictors) %*% lda.fit$scaling
matrix_product <- as.data.frame(matrix_product)

plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color=train$crime)
#tried this:plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color=km$cluster)

The dimensions of the model_predictors dataset are 404 and 13. The first 3D-plot is created using plotply-package with colors assigned according to “crime” (categorical variable created earlier). In the second 3D-plot the colors should be assigned according to the clusters found using kmeans-method. I had an error in this part so that only the first 3D-plot could be created and the plots wont be compared here.